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Abstract. Rana uenoi, a common brown frog inhabiting South Korea and the Japanese Tsushima Island, has a long history 
of taxonomic discussion. Although R. uenoi was described in 2014 as being distinct from R. dybowskii, the taxonomic sta- 
tus and distribution of the South Korean clade of R. uenoi has remained uncertain, particularly with respect to the related 
species R. dybowskii. Considering the species’ importance as a biological indicator of climate change, resolving the taxo- 
nomic equivocacy and clarifying its genetic structure is a matter of urgency. To do this, we included samples from across its 
distribution in South Korea, sequenced and analyzed Cytochrome b sequences, and developed novel microsatellite mark- 
ers for population genetic analyses. Our phylogenetic analyses verified the absence of R. dybowskii or another species in 
South Korea and confirmed the genetic divergence between R. uenoi and R. dybowskii. Population genetic analyses showed 
two distinct groups within R. uenoi, one on the Korean mainland and Japanese Tsushima Island and the second on Jeju 
Island, with the likely scenario being the Jeju Island population originating from mainland Korea during the Pleistocene. 
Demographically, we found evidence of population bottleneck events. Based on these results, we propose the English com- 
mon name “Korean large brown frog” and a Korean common name “241714 =|” [keun-san-gaeguri] for R. uenoi. These 
results provide an important baseline data for understanding future climate change impacts. 
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Introduction 


Protruding from Northeast Asia into the North Pacific 
Ocean, the Korean Peninsula functioned as a glacial refu- 
gium for numerous species during the Pleistocene glacial 
periods (CHUNG et al. 2018, LEE et al. 2018), especially for 
amphibians (ZHANG et al. 2008, BORZÉE et al. 2017, JEON et 
al. 2021). This function resulted in the peninsula harbour- 
ing today a significant biodiversity and cryptic endemism 
relative to its size (BAEK et al. 2011, CHUNG et al. 2018). In 
this respect, there has been a long-lasting discussion about 
the cryptic diversity and taxonomic status of a common 
brown frog species in South Korea, Rana uenoi, particu- 
larly in relation to R. dybowskii (MATSUI 2014, SONG & 
CHANG 2018). 


Since the establishment of binomial nomenclature and 
the description of the first brown Rana frog from Northeast 
Asia, the scientific name of the South Korean clade of R. uenoi 
has changed several times. This South Korean clade was first 
treated as a species of the Rana temporaria LINNAEUS, 1758 
group that was presumed to be distributed across the Eur- 
asian continent (STEJNEGER 1907). The clade was later as- 
signed to R. dybowskii GUNTHER, 1876 (GUNTHER 1876), fol- 
lowed by the its designation as a subspecies, R. temporaria 
ornativentris (OKADA 1928, MATSUI 2014) in 1928, and later 
changed to R. t. temporaria (SHANNON 1956). The taxonomic 
status of the clade in South Korea subsequently underwent 
a period of back-and-forth, oscillating between two subspe- 
cies of the Rana temporaria group: R. t. dybowskii, and R. t. 
chensinensis. The name R. dybowskii was used for the clade 
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in South Korea for the first time in 1978 (YANG & YU 1978) 
and was also confirmed as a species-level taxon in 1983 ow- 
ing to different chromosome numbers compared to R. tem- 
poraria (GREEN 1983, FROST 1985). 

Since then, R. dybowskii became the recognized name 
of the South Korean clade (Kim et al. 1999, 2002, LEE et al. 
2011) until Marsut (2014) assigned clades from the Japa- 
nese Tsushima Island and South Korea to the newly de- 
scribed R. uenoi. In doing so, Matsut (2014) split the new 
species from R. dybowskii from northern regions in Chi- 
na and Russia. Despite the description of the new species 
R. uenoi, this name has not yet gained general acceptance, 
and the denomination R. dybowskii is still frequently used 
in Korea (e.g., LEE & PARK 2016, Do et al. 2018, National In- 
stitute of Biological Resources 2019, The Korean Society of 
Herpetologists 2020, PARK et al. 2021), with some authors 
treating R. uenoi as a synonym of R. dybowskii (e.g., LEE & 
PARK 2016, National Institute of Biological Resources 2019, 
The Korean Society of Herpetologists 2020). The absence 
of a Korean common name for R. uenoi is another reason 
for this taxonomic equivocacy, and the Korean common 
name of R. dybowskii has been used for both species to 
date, despite the fat that there are even two Korean com- 
mon names for it: "YAFA T 2] [bukbang-san-gaeguri] 
and àZ) TE] [san-gaeguri] (SONG & CHANG 2018). Ad- 
ditionally, some herpetologists also suspect the presence of 
R. dybowskii and an undescribed cryptic species in South 
Korea due to the morphological variations within the 
South Korean clade of R. uenoi (LEE & PARK 2016). 

Rana uenoi is a common brown frog in South Korea, but 
its distribution is limited to the Korean Peninsula and near- 
by islands, including the South Korean Jeju Island and the 
Japanese Tsushima Island. As a result of overharvesting for 
food and the loss of breeding habitat, its population has been 
shrinking and it was consequently added to the list of spe- 
cies whose capture is prohibited in South Korea (National 
Institute of Biological Resources 2019). In addition, this spe- 
cies was listed as a biological indicator of climate change by 
the Korean government (National Institute of Biological Re- 
sources 2019), due to its breeding season shifting towards an 
earlier date, putatively due to global warming and increased 
winter temperatures (Kwon et al. 2020). In Korea, the ge- 
netic structure of R. uenoi populations has been studied us- 
ing mitochondrial and allozyme data, but these studies are 
limited in their geographic scopes (Kim et al. 1999, 2002). 
Further population genetic studies are necessary to establish 
a baseline and monitor changes before the species might be- 
come significantly impacted by future climate change (e.g., 
mortality of egg masses due to abnormal temperature fluc- 
tuations during warmer-than-normal winters). 

Here, we review the taxonomic status of R. uenoi using a 
DNA barcode gene Cytochrome b (Cyt b), define its popu- 
lation structure using Cyt b and nine newly developed mi- 
crosatellite markers, and clarify the presence or absence 
of R. dybowskii and any other related species in South Ko- 
rea. Contrary to previous studies using a limited number 
of samples (Kim et al. 1999, 2002, MATSUI 2014), here we 
analyzed data from samples collected across South Korea. 
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Materials and methods 
Sample collection 


Putative Rana uenoi samples were collected between 2006 
and 2020 from ten geographic populations in South Ko- 
rea. By including GenBank sequences from the type local- 
ity (Tsushima, Japan), we covered the whole distribution 
range of this species (Fig. 1 and Table 1). For tissues, we 
clipped a toe at the height of the last joint for adults or cap- 
tured tadpoles (from separated schools) and stored them 
in 70% EtOH until DNA extraction. Each adult was re- 
leased at the point of capture immediately after sample col- 
lection. Sampling was permitted by local governments ac- 
cording to the Wildlife Protection and Management Act of 
the Korean Ministry of Environment. The tissue sampling 
protocol conformed to the guidelines of Seoul National 
University Institutional Animal Care and Use Committee 
on ethical animal experimentation. We collected a total of 
239 samples for this study (see Table 1 for the detailed sam- 
ple sizes of each analysis). 


Laboratory protocols 


We extracted Genomic DNA with the DNeasy Blood and 
Tissue Kit (QIAGEN, Hilden, Germany) according to the 
manufacturer’s manual. Extracted DNA quantity and qual- 
ity were evaluated using an Epoch Microplate Spectropho- 
tometer (BioTek, Winooski, VT, USA). We uniformly di- 
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Figure 1. Sampling localities of Rana uenoi used in this study 
throughout South Korea. Copyright of the base map belongs to 
OpenStreetMap contributors and the GIS User Community. Yel- 
low dots indicate sampling sites for mitochondrial data, while 
red stars indicate sampling sites for microsatellite data. Green 
triangles indicate approximated localities for mitochondrial data 
of individuals from YANG et al. 2017. 
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Table 1. Sample sizes for mitochondrial DNA and microsatellite analyses of Rana uenoi populations in South Korea and the type 
locality (Tsushima, Japan; sequences from GenBank). 


Sample size 
Category Province City / County 
Phylogenetic Tree Haplotype Network Microsatellite Analyses 


14 


lon 
lon 


Mainland Gyeonggi Yangju 
Seoul 
Incheon 
Gangwon Chuncheon 
Yangyang 
Inje 
Gangneung 
Jeongseon 
North Chungcheong Goesan 
Cheongju 
Boeun 
South Chungcheong Asan 
Cheonan 
Seosan 
Yesan 
Gongju 
North Gyeongsang Pohang 
Gimcheon 
Daegu 
Gyeongju 
South Gyeongsang Hapcheon 
Hamyang 
Ulsan 
Yangsan 
Uiryeong 
Sancheong 
Haman 
Gimhae 
Jinju 
Busan 
Hadong 
Sacheon 
Goseoung 
Tongyeong 
Geoje 
North Jeolla Namwon 
Jeonju 
South Jeolla Jangseong 
Gurye 
Gokseoung 
Gwangyang 
Suncheon 
Jangheung 
Haenam 
Wando 
Island Jeju Jeju 


20 


Seogwipo 
Japan Tsushima Sago 
Shitaru 


Total 49 123 
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luted successfully extracted DNA samples to 10-20 ng/ul 
and stored at -20°C until the genetic analysis. 

We amplified the partial mitochondrial Cyt b region us- 
ing the primers Cytb 981F and Cytb 981R (Min et al. 2004). 
The PCR mixture consisted of total volume of 30 ul con- 
taining 1X Taq buffer with 2 mM MgCl, 0.2 mM dNTP 
mixture, 1 ul of 10 uM forward and reverse primers each, 
1 U Ex Taq polymerase (Takara Bio, Shiga, Japan), and 1 ul 
of template DNA using a TaKaRa PCR Thermal Cycler 
Dice Gradient. The PCR cycling procedure started with 
5 min at 95°C for initial denaturation; followed by 35 cycles 
of 40 sec at 95°C for denaturation, 40 sec at 45°C for primer 
annealing, 1 min at 72°C for extension; finishing with 7 min 
at 72°C for final extension. PCR products were sequenced 
on an AB 3730xl DNA analyzer (Applied Biosystems, Fos- 
ter City, CA, USA) by Cosmo Genetech Inc. (Seoul, South 
Korea). 

Novel microsatellite markers were developed for Rana 
uenoi. We outsourced whole genome paired-end sequenc- 
es of Rana uenoi (see SUK et al. 2021 for mitochondrial 
whole genome information), using the Illumina MiSeq 
platform (Illumina, San Diego, CA, USA), to Macrogen 
Inc. The adaptor sequences were removed from reads, then 
sequences were scanned to identify repeat DNA sequenc- 
es using RepeatModeler (Smit et al. 2014) (http://www.re- 
peatmasker.org/). We used SSR Finder (STIENEKE & Eu- 
JAYL 2019) to annotate reads containing simple sequence 
repeats. A total of 40 candidate markers (10 tetra-, 15 tri- 
and 15 dinucleotide repeats) were tested (amplification and 
polymorphism) on eight individuals (two samples from 
each of four geographic populations). Finally, a total of 13 
microsatellite markers were chosen to genotype all the re- 
maining samples. The microsatellite PCR mixture consist- 
ed of 1X Taq buffer with 2 mM MgCl, 0.2 mM dNTP mix- 
ture, 0.5 ul of 10 uM forward and reverse primers each, 1 U 
i-StarTaq polymerase (iNtRON Biotechnology, Seongnam, 
South Korea) and 1 ul of template DNA in a total volume 
of 20 pl. PCR cycling was conducted using a Takara PCR 
Thermal Cycler Dice Gradient (Takara Bio) with touch- 
down-PCR conditions consisting of an initial denaturing at 
94°C for 5 min, 20 cycles of denaturing at 94°C for 20 sec, 
touchdown annealing at 60-50°C for 20 sec, and an exten- 
sion at 72°C for 20 sec; followed by 20 cycles of denaturing 
at 94°C for 20 sec, annealing at 50°C for 20 sec, and an ex- 
tension at 72°C for 20 sec; followed by a final extension at 
72°C for 7 min. We outsourced the genotyping of amplified 
PCR products on an ABI 3730xl platform to NICEM Inc. 
(Seoul, South Korea). Microsatellite and mitochondrial se- 
quences used in this paper have been deposited at Gen- 
Bank (Accession Nos. MW 448298-448330, MW58891- 
58902, MW626913-626915). 


Cyt b phylogenetic tree and haplotype network 
We checked the quality and edited Cyt b sequence data 
using Geneious Prime (https://www.geneious.com). The 


sequences were aligned in MEGA-X (Kumar et al. 2018) 
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using the Clustal-W (THOMPSON et al. 1994) method un- 
der the default setting. We generated Cyt b haplotype data 
using DnaSP 5.10.01 (LIBRADO & Rozas 2009) and de- 
posited the haplotype sequences at GenBank (Accession 
Nos. MW448298-448330, MW58891-58902, MW626913- 
626915; Supplementary Table S1). 

For the phylogenetic tree reconstruction, we included 
Cyt b sequences of Rana frogs from GenBank, mostly Eur- 
asian species, including three R. uenoi samples from the 
type locality (27 species, 35 sequences), and used Pelophy- 
lax nigromaculatus from GenBank as an outgroup (Sup- 
plementary Table S1). The substitution model was selected 
according to the Bayesian Information Criterion computed 
by the greedy algorithm (LANFEAR et al. 2012) using Parti- 
tion Finder 2.1.1 (LANFEAR et al. 2017) before tree recon- 
struction. We reconstructed the Bayesian inference (BI) 
phylogenetic tree in MrBayes 3.2.6 (RONQUIST & HUELSEN- 
BECK 2003), treating missing values as non-contributing. 
The software conducted two independent runs of 107 Me- 
tropolis Coupled Markov Chain Monte Carlo (MCMC) 
generations using four chains (three heated and one cold 
chain; temperature set to 0.1) per run. Trees were sampled 
every 500 generations and summarized after discarding 
25% as burn-in. We visualized the final tree in FigTree 1.4.3 
(http://tree.bio.ed.ac.uk/software/figtree/). 

We reconstructed the haplotype network computed 
from the Cyt b sequences of R. uenoi using the median- 
joining approach in NETWORK 10.0.0.0 (BANDELT et al. 
1999) (http://www.fluxus-engineering.com). The network 
calculation relied on variable sites within the 731 bp of the 
Cyt b region from a total of 116 sequences, including three 
R. uenoi samples from Tsushima Island. 


Population genetic analyses using microsatellites 


We checked the quality of genotyping and peak calling us- 
ing Geneious Prime. We used Micro-Checker 2.2.3 (VAN 
OOSTERHOUT et al. 2004) to detect experimental errors, 
such as null alleles, large allelic dropouts, and scoring er- 
rors. We examined the likelihood of linkage disequilibrium 
between the 13 selected marker loci and the Hardy-Wein- 
berg equilibrium using Fisher’s exact test as incorporated 
in GENEPOP 4.7.2 (RoussET 2008) with default settings. 
We used GenAIEx 6.503 (PEAKALL & SMOUSE 2006) to cal- 
culate summary statistics for each locus and geographic 
population, such as number of alleles (N), number of ef- 
fective alleles (N_), observed heterozygosity (H,), expected 
heterozygosity (H,), and fixation index (F,.) for each lo- 
cus and population. Pairwise-F,,. (WRIGHT 1965) and -R,,. 
(SLATKIN 1995) values between populations were calcu- 
lated with 1,000 permutations using ARLEQUIN 3.5.2.2 
(EXCOFFIER et al. 2005). We used GenAIEx to carry out a 
Principal Coordinate Analysis (PCoA) with a covariance- 
standardized approach to identify genetic clusters among 
populations based on the distribution of individuals. We 
used STRUCTURE 2.3.4 (PRITCHARD et al. 2000) to in- 
fer the Bayesian population genetic structure from 10 in- 
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dependent runs of 10° MCMC generations (including 10% 
burn-in) from K = 1 to K = 10. The optimal K was deter- 
mined based on the EVANNo method (EvANNO et al. 2005) 
on STRUCTURE HARVESTER web 0.6.94 (EARL & VON- 
HOLDT 2012). Lastly, we searched for signs of bottleneck 
events using the mode-shift indicator and Wilcoxon signed 
rank test in BOTTLENECK 1.2.02 (CORNUET & LUIKART 
1997) with 1,000 iterations at TPM (0% SMM, 36% vari- 
ance; as recommended for most microsatellites by the soft- 
ware authors) and a Garza-Williamson index (M-ratio; 
GARZA & WILLIAMSON, 2001), which is lower than 0.68 
if populations have undergone population declines, using 
ARLEQUIN. 

For further analyses, we divided the Korean samples 
into two groups: one limited to Jeju Island (here referred 
to as “Island”), and one restricted to the Korean mainland 
and associated islands that used to be in direct contact 
with the mainland until a few thousand years ago (Li et 
al. 2014) (referred to as “Mainland”). The Japanese sam- 
ples were not used in this analysis. To understand the de- 
mographic history of Mainland and Island populations, 
two dispersal scenarios (Scenario 1 and Scenario 2) were 
compared using the approximate Bayesian computation 
incorporated in DIYABC 2.1.0 (CorNUET et al. 2014) with 
microsatellite data. Scenario 1 describes ‘from Mainland 
to Island’ dispersal and Scenario 2 ‘from Island to Main- 
land’ dispersal. The parameter settings for the two groups 
and the two scenarios were based on preliminary analyses 
that had followed the user manual. To test our scenarios, 
we applied a generalized stepwise mutation model with 


Pelophylax nigromaculatus 


0.08 


the generation time set to 2-3 years. Ranges of prior dis- 
tributions of parameters were adjusted from preliminary 
runs. 


Results 
Phylogenetic tree using Cytochrome b 


We recovered a 802-bp aligned segment including sites 
with missing characters for the phylogenetic tree analysis 
and a 731-bp aligned segment without missing characters 
for the haplotype network analysis. Among the 120 Rana 
uenoi individuals sampled in South Korea, 45 haplotypes 
were detected with DnaSP when including sites with miss- 
ing characters and 33 haplotypes when excluding sites with 
missing values. The best partitioning scheme deduced 
from these haplotypes was a 3-partition scheme according 
to codon positions: F81 (FELSENSTEIN 1981) + I for the first 
position, TRNEF (TAMURA & NEI 1993) + I + G for the sec- 
ond position, and TIM + G for the third position of co- 
dons. In the BI tree, all three R. uenoi individuals from Tsu- 
shima Island clustered with the Mainland populations of 
South Korea, having the same single haplotype as individu- 
als from Busan and Yangsan (in South Gyeongsang Prov- 
ince) which are the regions of the Korean Peninsula clos- 
est to Tsushima Island. There were no particular groupings 
within the Mainland populations in relation to geographic 
locations (Fig. 2). However, the Island populations in Jeju 
Island clustered separately from the Mainland populations 
(Fig. 2). 


Rana uenoi (Mainland populations, South Korea; Tsushima Island, Japan) 


Rana uenoi (Island populations, South Korea) 


R. dybowskii (Meilongjiang Province, Jilin Province, Shanxi Province, China; 
Primorsky Krai, Russia) 

R. chensinensis e 

R. kukunoris © s 

R. hecksheri 

R. zhengi 

R. arvalis 

R, asiatica 


- R. i 
R. longicrus B 
pirimi R. zhenhaiensis 
R. culaiensis i s 
p z R. dabieshanensis 
R. omeimontis p  muviensi 
R. hantuica PORET ae 
4 > R. chaochiaoensis 
R, tsushimensis 
R. uima 
R. amurensis 
C 
Node posterior probability 0.00 1.00 


R. grylio 
R. shuchinae 
R. johnsi 

R. temporaria 


R. coreana 


Figure 2. Bayesian phylogenetic tree of Rana uenoi and other species in the family Ranidae based on a fragment of the mitochondrial 
Cytochrome b gene. The Bayesian posterior probability (BPP) of each node is indicated by the gradient bar below the tree, where dark 
red represents a higher BPP. The green triangle indicates the position of the haplotype including Tsushima Island samples. Detailed 
information on populations of Rana uenoi and outgroup species can be found in Table 1 and Supplementary Table S1. 
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The haplotype network indicated a pattern concord- 
ant with that of the phylogenetic tree; The network did not 
highlight any particular grouping based on geographic lo- 
cation, but had one major haplotype including individuals 
from all eight administrative provinces of South Korea and 
all derivative haplotypes clustered around it (Fig. 3). Indi- 
viduals from Tsushima Island were characterized by the 
same haplotype as individuals from Busan, and northern 
Gyeongsang, southern Gyeongsang, and southern Jeolla 
Provinces. The cluster of haplotypes from Jeju Island in- 
dividuals was separated from Mainland haplotypes by 11 
mutational steps (Fig. 3). 
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Figure 3. The haplotype network of populations of Rana uenoi 
based on a fragment of the mitochondrial Cytochrome b gene. 
Circle sizes are proportional to the number of samples included 
in the haplotype analyses. Sites were grouped into provinces and 
each province was colour-coded. Detailed information on sample 
sizes from each locality, localities included in each province, and 
localities included in each haplotype analysis can be found in 
Table 1 and Supplementary Table S1. 
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Table 2. Descriptive indices of the nine microsatellite markers 
newly developed for Rana uenoi in South Korea. Na = the num- 
ber of alleles, Ne = the effective number of alleles, H, = observed 
heterozygosity, H, = expected heterozygosity, and F, = fixation 
index of each microsatellite marker. 


Locus Bae a Na Ne Ho H, Fis 
RuMic_4 ATAC 70-162 21 6.638 0.509 0.849 0.400 
RuMic_7 AGAT 172-484 49 27.604 0.761 0.964 0.211 
RuMic_10 AAAG 172-336 19 10.662 0.595 0.906 0.343 
RuMic_17 ATA 189-252 21 11.475 0.778 0.913 0.148 
RuMic_21 TG 104-250 43 24.669 0.558 0.959 0.419 
RuMic_28 GT 146-178 12 5.583 0.521 0.821 0.365 
RuMic_34 TAT 120-201 17 7.761 0.407 0.871 0.532 
RuMic_35 TAA 205-262 18 6.505 0.613 0.846 0.275 
RuMic_39 AT 101-143 12 4.417 0.571 0.774 0.262 


Population genetic analyses using microsatellites 


Thirteen microsatellite markers (7 tetra-, 3 tri-, and 3 dinu- 
cleotides) were developed out of the 40 candidates (depos- 
ited in GenBank, accession Nos. MW273296-MW273308; 
Supplementary Table S2). Among the 13 microsatellite 
markers, nine pairs exhibited evidence of linkage disequi- 
librium in GENEPOP (Supplementary Table S2), resulting 
in the exclusion of four markers (RuMic_2, 3, 6, and 9). The 
remaining nine markers (3 tetra-, 3 tri-, 3 dinucleotides) 
were used for further analyses. All markers had moder- 
ate levels of null alleles based on the results provided by 
MICRO-CHECKER, but we kept all nine markers because 
phylogenetically similar species are known to show similar 
characteristics (MARTINEZ-SOLANO et al. 2005, ZHANG et 
al. 2010). When we detected stuttering, peak calling was 
manually verified. Allelic dropout was not detected. 
Detailed information and descriptive indices for the 13 
microsatellite markers are summarized in Table 2 and Sup- 
plementary Table S2. The number of alleles ranged from 
12 (RuMic_28) to 49 (RuMic_7). Observed heterozygosity 
ranged from 0.407 (RuMic_34) to 0.778 (RuMic_17) and 
expected heterozygosity ranged from 0.774 (RuMic_39) to 
0.964 (RuMic_7). These discrepancies between H, and H, 
led to high levels of fixation indices. Descriptive indices of 
populations computed from these microsatellite markers 
are summarized in Table 3. The range of number of alleles 
spanned from 5.444 (Jeju) to 12.111 (Namwon). Observed 
heterozygosity ranged from 0.519 (Sancheong to 0.634 (JS) 
and expected heterozygosity ranged from 0.514 (Jeju) to 
0.859 (Boeun). The discrepancy between H, and H, for all 
populations, with the exception of Jeju and Seogwipo that 
had little discrepancy, led to high-level fixation indices. 
Pairwise-F,, and R,,. followed a similar pattern. Gener- 
ally, moderate to high levels of genetic differentiation were 
identified, and the populations on Jeju Island (Jeju and 
Seogwipo) consistently exhibited higher levels of genet- 
ic differentiation than any of the other Mainland popula- 
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Table 3. Descriptive indices of ten Rana uenoi populations in South Korea analyzed using microsatellites. Data shown for each popula- 
tion are the samples size (N), the number of alleles (Na), the effective number of alleles (Ne), observed heterozygosity (H,), expected 
heterozygosity (H,), unbiased expected heterozygosity (uH,), fixation index (F), mode-shift indicator, and Wilcoxon Test P (one- 
tailed probability for heterozygosity excess) from BOTTLENECK, and M-ratio. 


Popin N Na Ne Ho H, uH, Fs Mode-shift Wilcoxon M-ratio 

(SE) (SE) (SE) (SE) (SE) (SE) Indicator Test P (SE) 
Fangu 14 (1993) (0865) (0e) (0032) (00033) (oops) Tebaped 000293 (0178) 
Chuncheon 20 G39) (o)  @01) (0020 (002) (0.065) Shaped 010156 (0162) 
Parii 13 osa) (0623) (0055) (011) (0o) (oos) TSPePed 000098 A 
Usati 16 (1402) (253) (0.040) (0034) (0.035) (aoso EShaPed 012500 C0741) 
Sancheong 9 0716) (0938) (0069) (0030) (0032) (0o92) ‘Shaped 063281 (0.119) 
Geoje 1g (oes) (0.660) (0.081) (0.025) (one) (0.099) ee Ee (0.068) 
wmon 7 RUD BIA aaa oss aso O32 es gga BIS 
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Table 4. Pairwise-F „ and R,,.among 11 Rana uenoi populations in South Korea analyzed using nine microsatellites. Estimates of F, 
appear above the diagonal and estimates of R„„ below diagonal. All estimates significantly deviated from zero (p < 0.05) except those 
denoted by ‘NS: Abbreviations - YJ: Yangju, CC: Chuncheon, BE: Boeun, US: Ulsan, SC: Sancheong, GJ: Geoje, NW: Namwon, JS: 


Jangseong, JJ: Jeju, SGP: Seogwipo. 


YJ CC BE US SC GJ NW JS JJ SGP 
YJ 0.035 0.053 0.046 0.081 0.075 0.048 0.035 0.277 0.247 
CC 0.00555 0.037 0.045 0.048 0.062 0.031 0.038 0.262 0.233 
BE 0.04655 0.02755 0.054 0.057 0.056 0.01885 0.031 0.256 0.230 
US 0.01355 0.02755 0.006 NS 0.051 0.076 0.030 0.058 0.264 0.233 
SG 0.177 0.150 0.07155 0.133 0.081 0.02785 0.046 0.288 0.261 
GJ 0.271 0.252 0.191 0.267 0.122 0.042 0.031 0.260 0.238 
NW 0.232 0.188 0.118 0.183 -0.0335 0.072 0.021 0.275 0.241 
JS 0.153 0.121 0.0735 0.125 -0.0065 0.0335 -0.02055 0.247 0.221 
JJ 0.091 0.176 0.129 0.084 0.244 0.301 0.270 0.190 0.041 
SGP 0.117 0.219 0.177 0.117 0.336 0.423 0.369 0.283 0.01258 


tions. In the PCoA, Jeju and Seogwipo populations were 
clearly separated from all the other Mainland populations 
by principal coordinate 1 (Fig. 4). Similarly, the STRUC- 
TURE analysis identified two genetic clusters according 
to the delta K method (EvANNO et al. 2005), which distin- 
guished Island from Mainland populations (Fig. 5). 
Different signs of demographic history were indicated 
by the mode-shift indicator and Wilcoxon signed rank test 
and M-ratio (Table 3). The mode-shift indicator did not 
suggest any population size reduction in contrast to the M- 


ratio that indicated population size reduction for all popu- 
lations, as an M-ratio lower than 0.68 is considered to in- 
dicate a reduction in population size (GARZA & WILLIAM- 
SON 2001). On the other hand, the Wilcoxon signed rank 
test detected significant bottlenecks in four populations: 
Yangju, Boeun, Geoje, and Jangseong. 

From the result of DIYABC, Scenario 1, ‘from Main- 
land to Island’ dispersal, was supported slightly better than 
Scenario 2 by both the direct approach and the logistic re- 
gression results (Supplementary Tables S3, S4 and Fig. S1). 
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The divergence dates back to 19,700 generations (39,400- 
59,100 years, based on a generation time of 2-3 years). 


Discussion 


In this study, we review the molecular taxonomic status 
of Rana uenoi, a common brown frog in South Korea. We 
build on the data from previous studies (MATSUI 2014, 
Yuan et al. 2016, YANG et al. 2017) by including samples 
from across Korea, and verify the absence of R. dybowskii 
from South Korea. Through the development of new mi- 
crosatellite markers for R. uenoi, we are able to characterize 
the genetic diversity and population structure of this spe- 
cies in South Korea. 

Our sampling included R. dybowskii samples from Chi- 
na (Changbai Mountain in Jilin Province) and Russia (Pri- 
morsky Krai), regions proximate to the Korean Peninsu- 
la, and the samples of the South Korean/Japanese R. uenoi 
clade proved distinct from them. This is consistent with 
BorZEE et al. (2021) from which R. uenoi is known to be 
present as far north as Pyongyang in North Korea. Rana 
dybowskii and R. uenoi are estimated to have diverged in 
the Miocene due to the orogenesis of the Changbai Moun- 
tain Range approximately 11.01-8.55Ma (YUAN et al. 2016, 
YANG et al. 2017). Our results, are consistent with this hy- 
pothesis, as all samples from south of the Changbai Moun- 
tain Range in our analyses were referable to R. uenoi. In ad- 
dition, we did not find evidence for the presence of a spe- 
cies other than R. wenoi in Korea and R. dybowskii. The fact 
that all three R. uenoi individuals from Tsushima Island 
were grouped together within the Mainland populations 
strongly supports the status and distribution of R. uenoi as 
suggested by MATSUI (2014). 

The Korean Society of Herpetologists still lists R. uenoi 
and R. dybowskii under the same Korean common name, 
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Figure 4. The principal coordinate analysis plot based on micro- 
satellite F,. values for populations of Rana uenoi. Coordinates 1 
and 2 explain 13.75 and 4.93% of the genetic differentiation be- 
tween populations, respectively. Abbreviations - YJ: Yangju, CC: 
Chuncheon, BE: Boeun, US: Ulsan, SC: Sancheong, GJ: Geoje, 
NW: Namwon, JS: Jangseong, JJ: Jeju, SGP: Seogwipo. JJ and 
SGP refer to populations on Jeju Island, while all others are from 
the Korean mainland. Colour coding of localities follows their 
provinces colour in Figure 3 (e.g., Black for YJ in this figure and 
Gyeonggi Province in Figure 3). 
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and based on our results we suggest to establish a new Ko- 
rean common name for R. uenoi, replacing the current 
name = ¥J4t7}-" 2] [bukbang-san-gaeguri] where =p} 
[bukbang] means northern distribution (which is more rel- 
evant to R. dybowskii) and A7} T- =] [san-gaeguri] refers 
to the genus Rana, with a new name: FAF} E] [keun- 
san-gaeguri]; where = [keun] means ‘large and reflects 
the species’ relatively large body size compared to other 
Korean brown frog species. In addition, as Matsui (2014) 
did not define an English common name for R. uenoi, we 
suggest the use of “Korean large brown frog’, considering 
its primary distribution on the Korean Peninsula. 

The genetic diversity of the South Korean clade of 
R. uenoi was higher in Mainland populations than in Is- 
land populations, but the inbreeding coefficient showed 
a contrary pattern. This implies that geographic popula- 
tions on the Korean mainland have been altogether genet- 
ically more isolated and diverse than the populations on 
Jeju Island, perhaps reflecting complex geographic barri- 
ers on the Korean mainland (JEON et al. 2021). Regarding 
the species’ population genetic structure, it is worth not- 
ing that while individuals from the Mainland populations 
are grouped together with individuals from Tsushima Is- 
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Figure 5. The delta K graph (a) and bar plot (K = 2) from the 
STRUCTURE result (b) based on microsatellite data, showing 
genetic differentiation between populations of Rana uenoi. The 
delta K graph indicates that K = 2 is the optimal number of popu- 
lations. Abbreviations - YJ: Yangju, CC: Chuncheon, BE: Boeun, 
US: Ulsan, SC: Sancheong, GJ: Geoje, NW: Namwon, JS: Jang- 
seong, JJ: Jeju, SGP: Seogwipo. JJ and SGP refer to populations 
on Jeju Island, while all others are from the Korean mainland. 


Resolving the taxonomic equivocacy and the population genetic structure of Rana uenoi 


land, individuals from the Island population form a cluster 
that is differentiated from the Mainland/Tsushima. This re- 
sult is consistent throughout all analyses: pairwise-F.,/R,,, 
phylogenetic tree, haplotype network, PCoA, and STRUC- 
TURE result. The DIYABC analysis inferred that the dif- 
ferentiated Island population was the result of immigration 
from the mainland, which began before the Last Glacial 
Maximum when the mainland and the islands were con- 
nected, between 33,000 and 14,000 years ago (CLARK et 
al. 2009). Temporal changes of connectivity between the 
mainland and the islands are likely related to the pattern 
of sea level change in the Yellow Sea. When the continen- 
tal shelf was maximally emerged circa 21,000 years ago (Lī 
et al. 2014), it facilitated connectivity between amphibian 
populations (BoRZEE et al. 2020). However, Jeju Island 
was separated from the Korean Peninsula by the combined 
large palaeorivers, creating a boundary difficult to cross 
for small amphibian species (ANGELONE & HOLDEREGGER, 
2009; LE Lay et al. 2015), as is also seen in the Vistula Riv- 
er in Poland, marking the boundary between Hyla arborea 
and H. orientalis (STÖCK et al. 2012). 

Similar patterns of genetic isolation or subspeciation 
on Jeju Island have been reported for other animals, too, 
including for mammals (Siberian Roe Deer [Caproulus 
pygargus; KOH & RANDI 2001; KOH et al. 2013; PARK et al. 
2014; LEE et al. 2015, 2016], Lexmann’s Shrew [Sorex cae- 
cutiens; Kou et al. 2012], Striped Field Mouse [Apodemus 
agrarius; YOON et al. 2004; Ox et al. 2013] and Asian Lesser 
White-toothed Shrew [Crocidura shantungensis; LEE et al. 
2018]), and other amphibians (Jeju Salamander [Hynobi- 
us quelpaertensis; SUK et al. 2020] and Japanese Tree Frog 
[Dryophytes japonicus; JANG et al. 2011]). Faunal divergence 
between mainland Korea and Jeju Island has been attribut- 
ed to the formation of a temporary land bridge between the 
Korean Peninsula and Jeju Island during the Pleistocene 
glacial epochs (Yoon et al. 2004, Ox et al. 2013, PARK et al. 
2014, LEE et al. 2018), and it is hypothesized that individu- 
als of Rana uenoi dispersed from the Korean mainland to 
Jeju Island when a land bridge formed (Woo et al. 2013, LEE 
et al. 2018), became isolated when this bridge disappeared, 
and the Mainland and Island populations began to diverge 
genetically (JOHNSON et al. 2000, JANG et al. 2011, PARK et 
al. 2014) by random genetic drift and/or adaptation to the 
island’s environmental conditions (Jo et al. 2012). 

Signs of bottlenecks in all geographic populations were 
detected by the Wilcoxon sign ranked test and M-ratio, 
but not by the mode-shift indicator. This incongruence 
between bottleneck indices has also been noted in sever- 
al other empirical studies (HUNDERTMARK & VAN DAELE 
2010, SWATDIPONG et al. 2010, GONG et al. 2015) and could 
be attributed to different calculation methods and the re- 
sultant time scales (CORNUET & LUIKART 1997, LUIKART et 
al. 1998, GARZA & WILLIAMSON 2001). The M-ratio com- 
pares the number of alleles to the range of alleles and gen- 
erally resolves along a timescale of several hundred gen- 
erations (GARZA & WILLIAMSON 2001). On the other hand, 
the Wilcoxon signed rank test is a non-parametric test to 
check for a significant loss of alleles and it reflects two to 


four times the ‘effective population size (N) generations 
(LurkarT et al. 1998, VALENZUELA-QUINONEZ et al. 2014). 
Lastly, the mode-shift indicator examines allele frequen- 
cy distribution with the idea that bottlenecks reduce the 
occurrence of low-frequency alleles and increases that of 
intermediate- and high-frequency alleles (LUIKART et al. 
1998), which may occur over a few dozen generations (VA- 
LENZUELA-QUINONEZ et al. 2014). Considering that this 
species’ N, should be much higher than a few dozen in- 
dividuals (Supplementary Table S3), the mode-shift indi- 
cator reflected only the most recent demographic history. 
Overall, it seems more likely that R. uenoi populations, and 
especially the Mainland populations, suffered from popu- 
lation size reduction a much longer time ago than a few 
dozen generations and have since recovered. This is sup- 
ported by a star-shaped haplotype network, a common sig- 
nature of rapid population expansion after population size 
reduction (SLATKIN & HUDSON 1991, HEWITT 2004). The 
smaller population size of the Island than on the Mainland 
in DIYABC may be due to the founder effect or the differ- 
ence in the number of samples used for Mainland and Is- 
land (124 Mainland samples, 39 Island samples). 

The characteristics of the markers developed in this 
study could be affected by this bottleneck event. A wide 
range of alleles compared to the number of alleles, as de- 
tected by the M-ratio, was found in three markers (e.g. 
RuMic_z7, 10, and 21; in Table 2). The relatively low levels of 
observed heterozygosity (Table 2) may have resulted from 
genetic drift or inbreeding that have increased the number 
of homozygotes during population expansions after popu- 
lation size reductions (CORNUET & LUIKART 1997, KLEIN- 
HANS & WILLOWS-Munro 2019). The nine marker pairs in 
linkage disequilibrium also suggest the presence of previ- 
ous bottleneck events, since reduced effective population 
size can cause superficial linkages between unlinked loci 
(HILL 1981, LEDIG et al. 1999, FREELAND et al. 2011). Simi- 
lar genetic consequences were identified in the sister taxon, 
R. dybowskii (ZHANG et al. 2010), and in R. iberica, a spe- 
cies once restricted to a refugium and now endemic to the 
Iberian Peninsula (MARTiNEZ-SOLANO et al. 2005). 

We reviewed the taxonomic status of Rana uenoi using 
Cyt b sequences and microsatellites to investigate the pop- 
ulation genetic structure of this common and large brown 
frog in South Korea. From the phylogenetic analyses, we 
verified the absence of R. dybowskii and other species in 
South Korea. Previous studies (MATSUI 2014, YANG et al. 
2017, SUK et al. 2021) have identified R. uenoi and R. dy- 
bowskii as distinct species, which we confirm after using a 
larger dataset including samples from across South Korea. 
Based on our results, we recommend the full replacement 
of instances of “Rana dybowskii” with “Rana uenoi” in the 
species lists of the Korean government and academic so- 
cieties (National Institute of Biological Resources 2019, 
The Korean Society of Herpetologists 2020). We also sug- 
gest changing the Korean common name of R. uenoi to 
ZAH TE] [keun-san-gaeguri] and establishing “Korean 
large brown frog” as its English common name, both of 
which were previously lacking. From the population ge- 
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netic analyses, we conclude that Mainland and Island pop- 
ulations are differentiated, and the Mainland populations 
have experienced a population bottleneck in the Pleisto- 
cene. 

The establishment of future management strategies for 
R. uenoi should be based on the taxonomy and genetic re- 
sults of this study. As an indicator species of climate change 
designated by the Korean Ministry of Environment, we 
should monitor for future ecological and genetic changes, 
not only for the management of this species due to its high 
inbreeding level, but also for that of other wildlife on the 
Korean Peninsula. In this regard, our results provide im- 
portant baseline data for comparison with future studies. 
Similarly, commercial imports of R. dybowskii into Korea 
should now be outlawed, as it is not a native species of Ko- 
rea. Continued imports of R. dybowskii may result in the 
establishment of invasive populations that potentially have 
negative ecological impacts. 
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